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PENALTY METHOD WITH Pl/Pl EINITE ELEMENT APPROXIMATION FOR 
THE STOKES EQUATIONS UNDER SLIP BOUNDARY CONDITION 

TAKAHITO KASHIWABARA, ISSEI OIKAWA, AND GUANYU ZHOU 


Abstract. We consider the Pl/Pl or Plb/Pl finite element approximations to the Stokes equations 
in a bounded smooth domain subject to the slip boundary condition. A penalty method is applied 
to address the essential boundary condition u • n = g on dQ, which avoids a variational crime and 
simultaneously facilitates the numerical implementation. We give 0{/i^/^+e^/^+/i/e^/^)-error estimate 
for velocity and pressure in the energy norm, where h and e denote the discretization parameter and the 
penalty parameter, respectively. In the two-dimensional case, it is improved to 0(h + -h 
by applying reduced-order numerical integration to the penalty term. The theoretical results are 
confirmed by numerical experiments. 


1. Introduction 


In this paper, letting Q C (A^ = 2, 3) be a bounded smooth domain, we consider the Stokes 
equations subject to the non-homogeneous slip boundary condition as follows: 


( 1 . 1 ) 


where m : 12 —)■ and p : H —)■ K are the velocity and pressure of the fluid respectively, and ly > 0 

is a viscosity constant. Moreover, / represents the given body force, g the prescribed outgoing flow 
on the boundary T := dfl, and r the prescribed traction vector on T in the tangential direction, 
with a(u,p) = —pi + + iyu)'^) being the Cauchy stress tensor associated with the fluid. The 

outer unit normal to the boundary F is denoted by n. The first term in (l.l)i is added in order to 


u — vAu + 'S/p = 

f 

in 


div u = 

0 

in 

n, 

u ■ n = 

9 

on 

dn, 

n ® n)a{u,p)n = 

T 

on 

dn, 


ensure coercivity of the problem without taking into account rigid body movements. We impose the 


compatibility condition between (1.1)2 and (I.II3 which reads 

( 1 . 2 ) J^gdj = 0. 

The slip boundary condition (0)3-(1H1)4 (or its variant the Navier boundary condition) is now 
widely accepted as one of the standard boundary conditions for the Navier-Stokes equations. There 
are many applications of the slip boundary conditions to real flow problems; here we only mention the 
coating problem and boundary conditions of high Reynolds number flow m- For more details on 
the application side of slip-type boundary conditions, we refer to Stokes and Carey [27] and references 
therein; see also John |16j for generalization combined with leak-type boundary conditions. 


In the present paper, our motivation to consider problem (1.1) consists in dealing with some mathe¬ 


matical difficulties which are specific to its finite element approximation. As shown by Solonnikov and 
Scadilov in |26] (see also Beirao da Veiga |3| for a generalized non-homogeneous problem), proving the 
existence, uniqueness, and regularity of a solution of does not reveal essentially more difficulty 
compared with the case of Dirichlet boundary conditions. Then one is led to hope that its finite element 
approximation could also be treated analogously to the Dirichlet case. However, it is known that a 
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naive discretization of (1.1), especially when a smoothly curved domain is approximated by a polyhe¬ 
dral domain leads to a variational crime in which we no longer obtain convergence of approximate 
solutions. 

Let us describe this phenomenon assuming 5 = 0 and considering piecewise linear approximation of 
velocity. In view of the weak formulation of the continuous problem, see ( |2.5[ ) below, a natural choice 
of the space to approximate velocity would be (we adopt the notation of Section]^: 


(1.3) 


Vhn = {vh &Vh : Vh-nh = 0 on r^^}, Tt := dflh 


where denotes the outer unit normal associated to F^. Now suppose that N = 2 and that any two 
adjacent edges that constitute F/i are not parallel. Then, one readily sees that Vhn above reduces to 
14 = 14 n HQ{nh)^- As a result, the finite element solution computed using 14„ is nothing but the one 
satisfying the Dirichlet boundary condition, which completely fails to approximate the slip boundary 
condition. For N = 3 or quadratic approximation (or whatever else), we may well expect similar 
undesirable reduction in the degrees of freedom that should have been left for the velocity components, 
which accounts for the variational crime. 

One way to overcome the variational crime is to replace the constraint in ( |1.3[ ) by {vh ■ n) (P) = 0 for 
each boundary node P. This strategy was employed by Tabata and Suzuki [SD] where is a spherical 
shell; see also Tabata [2H1I21]. Extension of the idea to the quadratic approximation was proposed by 
Bansch and Deckelnick in [1] using some abstract transformation Gh : ll/i —11 introduced by Lenoir [22]. 
For n of general shape, the exact values for n{P) or n o Gh{P) may not be available. In this regard, 
some average of n/j’s near the boundary node P can be used as approximation of those unavailable 
values. This idea was numerically tested by Bansch and Hdhn in [2] for N = 3 and by Dione, Tibirna 
and Urquiza [S] for N = 2 (with penalty formulation), showing good convergence property. However, 
rigorous and systematic evaluation of those approximations is non-trivial and does not seem to be known 
in the literature. Moreover, implementation of constraints like (vh ■ n){P) = 0 in a real finite element 
code is also non-trivial and requires special computational techniques (see e.g. Gresho and Sani [T31 p. 
540] and [TJ Section 5]), which are not necessary to treat the Dirichlet boundary condition. 

In view of these situations, in the present paper we would like to investigate a finite element scheme 
to ([TtJ such that: 1) rigorous error analysis can be performed; 2) numerical implementation is as easy 
as for the Dirichlet case. With this aim we adopt a penalty approach proposed by Dione and Urquiza [l0| 


(see also i) which, in the continuous setting, replaces the Dirichlet condition ( |1.1| )3 by the Robin-type 
one involving a very small number (called the penalty parameter) e > 0, i.e., (j{u,p)n-n+ ^{u-n — g) = 0 
on F. At the weak formulation level, this amounts to removing the constraint v ■ n = 0 from the test 
function space and introducing a penalty term -{u-n — g,v ■ n)r in the weak form. Our scheme transfers 
this procedure to the discrete setting given on ^lh', see (3.1) below. Since the test function space for 


velocity is taken as the whole 14 involving no constraints, this scheme facilitates implementation, which 
serves purpose 2) mentioned above. It is indeed simple enough to be implemented by well-known finite 
element libraries such as FreeFem++ m and FEniCS [23] . as is presented in our numerical examples. 

Let us turn our attention to the error analysis. The first error estimate was given by Verfiirth m 
who derived in the energy norm for g = t = 0. The same author proposed the Lagrange 

multiplier approach in [321133] for r = 0. Later, Knobloch [TH] derived optimal error estimates (namely, 
0{h) for linear-type approximation and for quadratic-type approximation) for 5 = 0 and for 

various combinations of finite elements satisfying the LBB condition, assuming the existence of better 
approximation of n than Uh- The convergence (without rate) under minimal regularity assumptions was 
proved by the same author in m- A different proof of 0(/i^A j_estimate for the P2/PI element was 
given by [T] for 5 = r = 0, assuming that n o Gh is known. The technique using Gh was then exploited 
to study the penalty scheme in uni, again for the P2/P1 element and for 5 = r = 0. 

In the present paper, we study the penalty scheme for the Pl/Pl element combined with pressure 
stabilization and also for the PIb/PI element. Our method to establish the error estimate is quite 
different from those of the preceding works mentioned above. First, we address the non-homogeneous 
boundary conditions (l.I) which were not considered previously. Second, concerning the penalty 
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scheme, we directly compare (m,p) and (uh,Ph), whereas Dione and Urquiza [10] introduced a penalized 
problem in the continuous setting, dividing the error estimates into two stages. 

Third, we define our error (for velocity) to be ||{( — Uh\\H^{nh)y where u may be arbitrary smooth 
extension of u. This differs from [52] and [TUTU] in which the errors were defined as Iju — M/iUprqnnnh) 
and ||m — Ufi o Gh\\H^{Q)^ respectively. In view of practical computation, our choice of the error fits 
what is usually done in the numerical verification of convergence when flji 7^ fl. Compared with the 
method of Knobloch m who also employed ||u — u/illprqnfc) as the error, the difference lies in the way 
the boundary element face S on T/i is mapped to a part of T. In fact, he exploited the orthogonal 
projection from T to T/i, the image of which is localized to each S. Then he needed delicate arguments 
(see m pp. 142-143]) to take into account the fact that it is not globally injective when N = 3 (this 
point seems to be overlooked in [52] p. 709]). We, to the contrary, rely on the orthogonal projection tt 
from T/i to r, which is globally bijective regardless of the space dimension N, provided the mesh size h 
is sufficiently small. This enables us to transfer the triangulation of to that on T in a natural way, 
which is convenient to estimate surface integrals. Complete proofs of the facts regarding tt used in this 
paper, which we could not find in the literature, are provided in Appendices P] and [B] 

Finally, we comment on the rate of convergence we obtain in our main result (Theorem |4.l[ ) 

which is not optimal. In our opinion, all the error estimates reported in the preceding works, which use 
Uh to approximate n, remain 0{h^^^). Verfiirth [321 Theorem 5.1] claimed 0{h); however, the estimate 




{uh-noTT ^)ah 


< C/i]]ull^i/2(r,^) 11 cTft, 1177-i/2(r^), 


which was used to derive equation (5.12) there, seems non-trivial because nh is not smooth enough 
globally on T^^ (e.g. it does not belong to i7^/^(r^)). If ||cr;i||7/-i/2(r^) on the right-hand side is replaced 
by l|f^ii||L2(rh)i then one ends up with 0(/i^/^) in the final estimate. Dione and Urquiza [M Theorem 
4] claimed 0(/i^/^); however, in equation (4.13) there, they did not consider the contribution 


1 




ll(«. 


Vh) ■ n\\L2(T), 


which should appear inside the infimum over Vh G Kh even when ^ (see Proposition 4.2 of Layton 
m)- If this contribution is taken into account, one obtains for the final result. 

To overcome the sub-optimality, in Section we investigate the penalty scheme in which reduced- 
order numerical integration is applied to the penalty term. This method was proposed in [5] and was 
shown to be efficient by numerical experiments for N = 2. We give a rigorous justification for this 
observation in the sense that the error estimate improves to 0{h) if e = 0{h?'). Our numerical example 
shows that the reduced-order numerical integration gives better results also for N = 3, although this is 
not proved rigorously. 

In our numerical results presented in Section [^ we not only provide numerical verification of con¬ 
vergence but also discuss how the penalty parameter e affects the performance of linear solvers. We 
find that too small e can lead to non-convergence of iterative methods such as GMRES, whereas sparse 
direct solvers such as UMFPACK always manage to solve the linear system. 


2. Formulation of the Stokes problem with slip boundary condition 

We present our notation for the function spaces and bilinear forms that we employ in this paper. 
The boundary F of D is supposed to be at least C^d.gniooth. The standard Lebesgue and Sobolev(- 
Slobodetskii) spaces are denoted by LP(D) and 1U®’P(D) respectively, for p G [l,oo] and s > 0. When 
p = 2, we let i7®(r2) := 1U'*’^(D). Their vectorial versions are indicated as LP(D)^ and etc.; 

however, in case they appear as subscripts, say jj • 11 lp(q)n, we write jj • llLp(n) for the sake of simplicity. 
The spaces above may also be defined for F; see e.g. [21]. 

We define function spaces to describe velocity and pressure as follows: 

U = Ffi(D)", Q = L2(D). 
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Also we set 

Vn ■= {u G : (Trz;) • n = 0 on F}, 

Q ■= = {qG L^{^) : f^qdx = 0 }, 

where Tr stands for the trace operator; for simplicity, Tr v is indicated as v in the following. Finally, to 
describe a quantity which corresponds to the normal component of the traction vector, i.e. a(u,p)n ■ n, 
we introduce 

( 2 . 1 ) = 

where the prime means the dual space. 

Let G C be an open set. The L^(G)- and L^(dG)-mner products are denoted by (•, •)g and (•, OoG, 
respectively. Moreover, we define bilinear forms ao, bo, cgc, for u,v £ q S L‘^{G), A,/r € 


L\dG), by 


(2.2) 

aciu^v) 

(2.3) 

bG{v,q) 

(2.4) 

CdG{X,p) 


- / (Vm + (Vuf) : (Vv + (Vvf) dx, 


(2.5) 


JdG 

where the dot and colon mean the inner products for vectors and matrices, respectively. When G = 
we use the abbreviation a = an, b = 6a, c = cga. In the following, c(', •) is also interpreted as the 
duality pairing between iJ^/^(r) a nd 

The variational formulation for (1.1) consists in finding (u,p) G V x Q such that u ■ n = g onT and 
I a{u,v) + b{v,p) = {f,v)n + ('r,r;)r Vn € 14, 

I 6(m, q) =0 yq G Q. 

It is well known that this variational equation admits a unique solution under the compatibility condition 
( |1.2[ ) and that it becomes regular according to the smoothness of T, /, g, r; see [3 US]. Throughout 
this paper, we assume T G G^’^, / G L^(n)^, g G iF^/^(r), r G , so that the regularity 

{u,p)G Hh^rxHHn) is assumed. 

Letting A := —a{u,p)n ■ n, we see that {u,p, X) G V x Q x A satisfies 

a(u, v) + 6(n, p) + c(n • n. A) = (/, v)n + (r, n)r Vn G V, 
b{u, q) =0 yq G Q, 

c{u ■ n — g, p) = 0 ypL G A. 


( 2 . 6 ) 


In fact, (2.6)i follows from Green’s formula. By (1.2) one has 6 (m, 1) = 0, which implies (2.6)2- Multi¬ 
plying u ■ n = g hy a test function p and integrating over F lead to (2.6)3. 


Remark 2.1. Observe that {u,p + k, X + k) with any A; G M is also a solution of (2.6). According to this 
fact, we will adjust the additive constant of p (and thus of A) later on, before we start error analysis of 
the finite element approximation (see Remark 4.1 below). 


3. Finite element approximation 


3.1. Triangulation and FE spaces. Let us introduce a regular family of triangulations Th oi a poly¬ 
hedral domain O/i, which is assigned the mesh size h > 0. Namely, we assume that: 

(HI) each T G Th is a closed A^-simplex such that Ht ■= diamT < h; 

(H2) rth = Utgt. T; 

(H3) the intersection of any two distinct elements is empty or consists of their common face of dimension 
< A^- 1; 

(H4) there exists a constant G > 0, independent of h, such that px < Ghx for ad T G Th where px 
denotes the diameter of the inscribed ball of T. 
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We define the boundary mesh Sh inherited from Th by 


Sfi = {S £ Sh ■ S is an {N — l)-face of some T e Th}- 


Then we see that Sh satisfies the requirements that are analogous to (H1)-(H4) above, and especially 
we have T/j := dflh = UseSh We assume that T/j approximates T in the following sense: 

(H5) the vertices of every S G Sh lie on T. 

Throughout this paper, we confine ourselves to the case where /i > 0 is sufficiently small, which will not 
be emphasized in the following. In particular, all the results given in Appendices [A] and [B] are supposed 
to hold true. 

As mentioned in Section we focus on the Pl/Pl and Plb/Pl finite element approximations for 
velocity and pressure, to which we refer as I = 1 and I = 16, respectively. Namely, we define 


Vh 


Qh 


I {vh e C{nh)^ : e Pi(t)^ vt e %} if i = h 

\{vh e c{nh)^ ■-vhW ^ Pi{T)^ ® b{t)^ yreTh] if i = 16 , 

{qh G CiUh) : qh\T& Pi{T) VT € %}, 


where B{T) stands for the space spanned by the bubble function on T. We also set Vh ■= 
and Qh ■= Qh H L'^{^h)- We consider a discrete version of the space A (recall (2.11) as 

kh = {\h&L\Th) : \h\s&Pi{S) V 5 G 54 , 


which is a discontinuous Pi finite element space on P/j. 

We turn our attention to interpolation operators. In order to deal with the situation T "we 
first extend functions defined in to which is a fixed bounded smooth domain containing Q U Qh- 
For this purpose we consider an arbitrary extension operator P : —)■ satisfying the 

stability condition 

||P/IU^,P(0) < C\\f\\w-,.^n) V/ G 

where C is a constant depending only on iV, 12, to, p- Such P does exist; for example, if to = 0 we may 
exploit the zero-extension as P, and if to > 1 it can be constructed e.g. by Nikolskii’s method (see ^^)- 
In view of the lift theorem which concerns a right inverse of the trace operator, we may also extend 
functions given on P to those in 12. We agree to use the same symbol P to refer to this extension as 
well. Then, given a function / in 12 or on P, we denote Pf by / for simplicity in the notation. 

Now we let Ih and Rh represent the Lagrange interpolation operator to linear PE spaces and a local 
regularization operator to linear PE spaces, respectively. Then, from the theory of interpolation error 
estimates (see [H Section 4]) combined with the stability of extension operators, we have 

II/ - 4/||H^(n,) < Ch^-n\f\\HHn), f G H^Q), to = 0,1, 

11/- i?J||H^(n,) < Ch^-n\f\\mm, f G TO = 0,1, 


where the constants C depend only on N, on the constant in assumption (H4) above, and on a reference 
element. Eor {Th)-estiuiates, we first apply a trace inequality on each boundary element and then 
add them up to obtain 

ll/-4/l|L^(r,) < C\\f-IhfYLi%Jf-hf\\]i^n.) 

< Ch^^^WfWHHT ^ Ch^^^WfWH^iD f G 


3.2. FE 

choose e 


scheme with penalty. We propose the following discrete problem to approximate (1.1): 
> 0 suitably small and find {uh,Ph) & Vh x Qh such that, for all {vh,qh) & Vh x Qh, 


ahiuh, Vh) + bh{vh,Ph) + -Ch{uh ■ Uh - hg, Vh ■ nh) 

e 

bh{uh,qh) 


{f,Vh)n^ + {f,Vh)r,,, 
dh{Ph,qh), 


(3.1) 
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where nu is the outer unit normal of and at ■= bh := bci^, Ch ■= cr^ (recall (2.2)-(2.4) above). 
Moreover, /, g, f represent the extensions of /, g, r respectively, which are discussed in the previous 
subsection, dhi-,-) is a pressure-stabilizing term, which is present only when I = 1 and is defined by 

1 if ; = 1, 


dh{Ph,qh) = rih^{Vph,^qh)ah 


T] := 


0 if I = lb. 


We remark that ry for the case I = 1 can be any positive constant; here, we suppose it to be 1 for 
simplicity. We state the well-posedness of this discrete problem. 


Proposition 3.1. There exists a unique solution (uh^Ph) G 14 x Qh of (3.1). 

The proof relies on the following discrete versions of Korn’s inequality and the inf-sup condition: 

(3.2) ot\\vh\\]ji(Q,^)<ah(vh,Vh) VvhGVh, 

(3.3) C\\qh\\LHnh)< sup nO,'”'"'*' -b Cgh\\\/qh\\LHQH) ^qh&Qh, 


C\\qh\\LHnf,)< sup + Cgh\\Wqh\\L^Un) 


where a > 0 depends only on TV, 12, u and C > 0 depends only on N^Tl. Proofs of these uniform 
coercivity estimates are found in na Theorem 4.3] and in [5S1 Proposition 4], respectively. They will 
be of central importance when we perform error analysis in Section]^ Now we prove the proposition. 


Proof of Proposition \3.1\ We notice that (3.1) is equivalently rewritten as follows: find {uh,Ph) G 14 x 
Qh such that, for all (vh, qh) & Vh x Qh, 


(3.4) 


Bh{uh,Ph',Vh,qh) ■=ah{uh,Vh) + bh{vh,Ph) - bh{uh,qh) + dh{ph,qh) + ^Ch{uh ■ nh,Vh ■ Uh) 

€ 

= if, Vh)nh + + -Chihg, Vh ■ nh). 


Because the problem is finite dimensional, it suffices to prove that Bh{uh,Ph',Vh,qh) = 0 for all Vh and 
qh implies Uh = Ph = 0. Taking {vh,qh) = {uh,Ph) yields, thanks to ( [T^ , Uh = 0 and dh{ph,Ph) = 0. 
Below we deal with I = 1 and I = lb separately. 

Let I = 1. Then, since 77 > 0 we have Vp/i = 0, which implies that ph equals some constant k. 
By (3.4) one has bh{vh,k) = —k Vh ■ Uhd'jh = 0 for all Vh G 14- Choosing Vh G 14 such that 
/p^ Vh ■ Uh d'jh 4 0 gives A: = 0, so that ph = 0. 

When I = lb, (3.4) reduces to b{vh,qh) = 0 for all Vh G 14. This combined with (3.3) tells us that 
Ph is equal to a constant. Discussing as in 1 = 1, we conclude ph = 0. This completes the proof of 
Proposition |3.1[ □ 


Let us rew rite (3.1) in such a way that the discrete problem becomes comparable with the continuous 
one given in (2.6). To this end, we introduce Xh '= \{uh ■ rih — Ihg) G A^, to obtain 


(3.5) 


ah{uh,Vh) + bh{vh,Ph) +Ch{vh ■ nh,Xh) = {f,Vh)n^ + {f,Vh)rh 

bh{uh,qh) = dh{ph,qh) 

Ch{uh ■ nh - Ihg, Ph) = ech{Xh, Ph) 


^Vh G Vh, 
yqh G Qh, 
yph G Ah. 


In fact, (3.5)i results from the symmetry Ch{Xh,Ph) = Ch{ph,Xh), and (3.5)3 follows from multiplying 
the equation Uh • Uh — Ihg = cA/i by any test function in L'^fTh) D A^. 


3.3. Auxiliary lemmas. We collect several results which will be useful to evaluate the difference of ft 
and Uh in terms of volume or surface integrals. We denote the symmetric difference of D and Qh by 

Q/xQh ■— {Q \ Qh) u {Qh \ Q) 

and call it the boundary skin. The first lemma concerns a linear operator which extends functions in I 4 
to Q is equipped with suitable stability properties. For the proof, we refer to na Theorem 4.1]. 
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Lemma 3.1. There exists a linear operator Ph -Vh - 
\\PhVh\\Hm(^a,) C\\vh\\H’^{nh) 

\\PhVh\\H^{QAQf,) < Ch^^‘^\\vh\\H^{nh) 
where the constants C are independent of h. 


such that PhVh\fi^ = Vh und 
'ivh eVh,m = 0,1, 
yvh eVh,m = 0,1, 


The second lemma gives estimates of volume integrals over the boundary skin, which are not restricted 
to discrete spaces. However, notice that we require higher regularity for the right-hand side. 

Lemma 3.2. Let / G to = 0,1. Then we have 

WfWH^iQAflh) ^ C'^ll/llH-+i(n)’ 

where C is independent of h. 

Theorem 


Proof. Since LlALlh C r(C'o£;h^) by Proposition A.2 Theorem A.3 for di = Cosh^ and for p = 2, 
combined with the trace inequality ||/||L 2 (r) < C'H/llHqn), leads to the desired result. □ 

The last lemma restates results concerning surface integrals obtained in Theorems |A.1| and |A.2[ 

Lemma 3.3. Let tt be the orthogonal projection to T defined in a tubular neighborhood o/T. Then we 
have the following estimates: 

||/o7r|U2(r,)<C||/|U2(r) V/S L2(r), 


fd^- fo'^d'jh 


< Ch^ 


iLi(r) 


v/eL'(r), 

yfGH\n), 


where the constants C are independent of h. Here, d'j and d'jh denote the surface elements associated 
withT andTfi, respectively. 

Remark 3.1. By the first and third estimates together with trace inequalities, for all / S 
(hence its extension / is in W{Ll)) we obtain 

||/||L2(rh) < II/- f °A\L'^(Th) + Il/0 7r||i2(r,^) < C||/||^i(5=2) < C\\f\\H^/2(j^). 

4. Error analysis of penalty FE scheme 

4.1. Estimation of consistency error. To explain the idea, suppose that Ll = Llh and that g = 0. 
Then subtracting the discrete problem (3.5) from the continuous one (2.6) would give us 

{ a{u - Uh, Vh) + b{vh,p - Ph) + c{vh ■ n, X - Xh) = 0 \/vh G Vh, 

b{u- Uh,qh) =-dh{ph,qh) ^qh G Qh, 
c{{u - Uh) ■ n, pLh) = -€c{Xh, fih) yph G Ah. 

This type of relation is essentially important in error analysis of the finite element method and is 
sometimes called the “Galerkin orthogonality”. If SI yf then such a relation is by no means available 
because subtraction is impossible. However, even in this case one can still expect that an asymptotic 
version of (4.1) should hold as h becomes small. The next proposition verihes this expectation. 

Proposition 4.1. Let {u,p,X) and (uh,Ph,Xh) be solutions of and l[3.5\ ) respectively. We assume 
the regularity of the data: T G C^’^, / G L'^(Ll)^, g G W/^{V), r G HP'^iT)^. Then there exist 
constants C = C{N,Q.,v, f,g,T), independent of h and e, such that for all {vh,qh, Ph) G Vh x Qh x Ah 
there holds 

\ah{u - Uh,Vh) + bh{vh,p - Ph) + Ch{vh ■n,X- A/i)| < Ch\\vh\\m(Q.,^), 

(4.2) { \bh{u-Uh,qh) + dh{ph,qh)\<Ch\\qh\\L‘2{nH)^ 

\ch{{u - Uh) ■ nh,Ph) + eCh{Xh,Ph)\ < Ch\\ph\\L^rh)- 
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Proof, (i) Let us prove (gi. Since “ /n\n,- obtain 

ah(u - Uh,Vh) = a{u,PhVh) - ah{uh,Vh) + aQ^\Q{u,Vh) - aQ\fi^{u, PtVh), 
bh{vh,p-Ph) = b{PhVh,p) - bh{vh,Ph) + bn^\n{vh,p) - bn\Q^{PhVh,p). 

On the other hand, it is clear that 

Ch{vh Xh) = c{PhVh • n, A) - Ch(vh ■ rih, Xh) + Ch{vh ■ rih, X) - c{PhVh ■ n, A). 


Addition of the three equations above combined with (2.6li and (3.5)i yields 

ah{u - Uh,Vh) + bh{vh,p-Ph) + Chivh ■n,X- Xh) 

= (/, PhVh)n - (/, Vh)Q^ + (r, PhVh)r - (f, Vh)rh 

+ anh\n{u,Vh) - PhVh) + bQ^\Q{vh,p) - bn\Q^{PhVh,p) 

+ Ch{vh ■ rih, A) - c{PhVh ■ n, A) 

(4.3) =:/i +/2 + /3 + /4 + /5. 


First we estimate volume integrals. It follows from Lemmas 3.1 and |3.2| together with the stability 
of the extensions that 

|4i| < \\f\\L^(nAnh)\\PhVh\\L^{nAn^) < Ch\\f\\L2(^Q)\\vh\\m(nh)’ 

|.?3| < C\\u\\Hi-{QAnh)\\Phyh\\H^{QAnh) ^ 

\h\ < C\\PhVh\\H^{nAnH)\\p\\L^{^i^iih) ^ Cb.’^^^Wplln 
Next we estimate surface integrals. For I 2 we observe that 

l2= T ■ PhVh dj- (t • PhVh) OTTd'Jh + {TOTr-f, {PhVh) o 7r)r,, 

Jr Jru 

+ (f, (PhVh) OTT- PhVh)rh 

=: I 21 + I 22 + l23- 


^53 


l54- 


It follows from Lemma 3^ Remark |3.1[ and Lemma |3.1| that 

I/ 21 I < Ch'^Wr ■ PhVh\\L^{r) < C'^^||u||L2(r) \\vh\\H^{a,^)i 
I/ 22 I < Ch\\f\\jji^^-^\\PhVh\\Hi(n) < Ch\\r\\Hi/2(r)\\vh\\H^nh)^ 

Ihsl < ||u||L2(r^) • Ch\\PhVh\\Hi(^Q-^ < Ch\\T\\fji/2(^r)\\r’h\\H^nH 
For /s we observe that 

I 5 = Ch{vh ■ {uh-no tt), A) + Ch{{vh - PhVh) ■ {n o tt), A) 

+ ChiiPhVh • n) o TT, A - A o tt) 

+ / {PhVh ■ n X) o TT djh - / PhVh ■ nXd'y =: I51 + I5 
Jvh Jr 

It follows from Lemmas |B.1[ |3.3[ and |3.1| that 

I/ 51 I < C'/i||?;;i||L 2 (r,^)||A||L 2 (r,^) < C'/i(||u||/i- 2 (o) + |b|l_H-i(n))lkh||/i-i(nh), 

\I52\ < \\PhVh - {PhVh) o7r||i2(r,^)||A||L2(r^) < Ch\\PhVh\\Hi(Q)\\M\L2{rh) 

< C'/i(||u||//2(n) + |b||H-i(n))lk/t||//i(ah)> 

l^sl < \\{PhVh ■ n) o 7r|ji2(r,^)||A - A o 7r||i2(r,^) < CWPhVhWL^^^r) ■ C'^llA|1 ^i(q) 

< Ch{\\u\\H2(Q) + lbllH-i(n))||u;j|j//i(a,^), 

b54| < Ch^\\PhVh ■ nAlbqr) < Ch'^\\PhVh\\L'^{r)\\X\\L^(T) 

< C'/i^(||u|l/i-2(Q) + |b||_Hi(n))|b/i||_Hi(nh)- 
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Combining the above estimates with (4.3) and noting that ||m||_h' 2 (q) + ||p||//i(n) < C(||/||L 2 (n) + 
llffllH3/2(r) + ||T||jji/ 2 (r)) by the regularity theory of the Stokes equations, we conclude (4.2)i. 


(ii) Let us prove ( | 4 . 2 [ ) 2 - One finds from ( | 2 . 6| )2 and ( 3 . 5)2 that 

bh(u - Uh,qh) + dh{ph,qh) = bh{u,qh) - b(u, Phqh) =■ h- 

By Lemmas |3.1| and |3.2| we have 

I/el < C\\u\\Hi{nAnh)\\Phqh\\L^(nAQH) < Ch^^^\\u\\H^(^n)\\qh\\L2(n^), 
from which ( 4 . 2)2 follows. 


(iii) Let us prove ( 4 . 2 ) 3 . One finds from ( 2 . 6)3 and u - n = g onT, which is due to ( 3 . 5 ) 3 , that 
Chiiu - Uh) ■ n/i, Hh) + echiXh, Ph) = Ch{u ■ Uh - hg, Ph) 

= Ch{u- [rih-no ■n),Hh) + Ch{{u - uo tt) ■ n o tt, ph) + Ch{g on- g, gn) 

+ Ch{g - Ih9,Ph) 

='■ -f? + /s + ^9 + -fio- 

By Lemmas B.l 3.31 and |3.1[ together with the interpolation error estimate, we have 
I/ 7 I < Ch\\u\\L-^(Tf^)\\ph\\L^(Th) < Ch\\u\\Hi{n)\\9h\\L-^{Th)^ 
l/sl < \\u - UO n\\L2(Y^)\\g,h\\L-^{Tn) < Ch\\u\\H^(n)\\gh\\L‘^{T,,), 

I/ 9 I < C'/i||5||^i(Q)||/r,i||i2(r,^) < C'/i||5||^i/2(r)||M/illL2(rh), 

|/io| < C'/i^/^||g||/i-2(Q^)||^,,||i2(r^) < C'/i^/^|jg||^3/2(r)||Ai?i||L2(rh), 

from which ( 4 . 2)3 follows. This completes the proof of Proposition 14. ll 


□ 


4.2. Error estimate for velocity and pressure. We are ready to state the main result of this paper. 

Theorem 4.1. Assume that L S C^’^, f € g S //^'^^(L), r G and that h > 0 is 

sufficiently small. We let {u,p) and {uh,Ph) be solutions of (2.5^ and l[3.1\ ) respectively. Then there 
exists a constant C = C(iV, O, u,/, ( 7 , r), independent of h and e, such that 

Wu-UhWHHn^) + \\{p + kh) - PhWL^n^) < C{Vh + ffe+ ^), 

where kh = - RhP, l)n;,. 

Remark 4.1. As mentioned in Remark |2.1[ there is room for us to choose arbitrary additive constant 
k for p. Here, we adjust k in such a way that Rh{p+ k) —ph G L'^fflh), which gives rise to the constant 
kh above. By considering {p + kh, A + kh) instead of (p. A), we may assume = 0 in the following proof. 


Proof of Theorem \4.1\ Let Vh '■= ffu, qh ■= RhP, Th '■= RhX. By interpolation error estimates one has 

l|w — Uh\\H^{nh) ^ C^\W\\H^(n)k + \\vh — u/tUnqnfe), 

\\p - PhWmn^) < C'lbllffba)^ + hh - PhWmn^), 

hence it suffices to bound \\vh - UhWmiQh) and \\qh -Ph\\m(nh) by C(h}^'^ + According 

to the uniform ellipticity (3.2) we obtain 

0 i\\vh - Mhll^qn,^) + dhijph - qh,Ph - qh) + £Ch{Xh - Ph,Xh - Ph) 

< ah{vh - Uh,Vh - Uh) + dh{ph - qh,Ph - qh) + eCh{Xh - Ph,Xh - Ph) 

= ahivh -u,Vh- Uh) 


+ ah{u - Uh, Vh - Uh) + bh{vh -Uh,p- Ph) + Ch{{vh - Uh) ■nh,X- Xh) 

- bh{vh - Uh,p-Ph) + dhijPh - qh,Ph - qh) 

- Ch{{vh - Uh) ■nh,X- Xh) + ech{Xh - Ph,Xh - Ph) 

=■ Ii + I 2 + I 3 + Ii. 


(4.4) 
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From interpolation error estimates and Proposition | 4 . 1 | it follows that 


(4.5) |/i| + I/ 2 I < gW'^h — + Ch^. 

For I 3 , we observe that 

I 3 = bh{vh -u + u- Uh,Ph - qh + qh-p) + dh{ph,Ph - qh) - dh{qh,Ph - qh) 
= bh{vh -u,qh-p) + bh{vh - u,ph - qn) + bhiu - Uh,qh - p) 

+ bhiu - Uh,Ph - qh) + dhiPh,Ph - qh) 

— dhiqh,Ph — qh) =■ I31 + I32 + I33 + I34 + l 35 - 

The first three terms are estimated as 

\l 3 l\<Ch^, \l 32 \<Ch\\ph-qh\\mn,), \l 33 \<Ch^ + ^\\vh-Uhrm^a,), 


whereas we know that IJ34I < Ch\\ph — qh\\L^(Qh.) Proposition 4.1 I35 is bounded, thanks to Holder’s 
inequality, by 

\l35\ < vb'^\Nqh\\h(^a^) + ^dhiph - qh,Ph - qh) < Ch^ + ^dhiph - qh,Ph - qh)- 

To obtain further estimates of I32 and I34, we observe from (| 3 . 3 |) that 


C\\ph-qh\\LH<:iH)< sup —^ + C'??/i||V(p,, 

. bhivh,Ph-p) , bhivh,p-qh) , ^ ,|| 

- “n—n-+ ®’^P “ 11 —li- +Cvh\\Wiph-qh)\\mn^)- 

VH&Vh \\'^h\\m{Q^) W'VhWmiQ,,) 

Here, the second term on the right-hand side is bounded by Ch\\p\\Hi (Q) = Ch. We claim that the first 
term is bounded by Ch + C\\vh — UhWn'^iQh)- In fact, it follows from ( 4 . 2 q, in which Vh is restricted to 
Vh (hence = 0 on P^t so that Chivh ■ Uh, ■) = 0), that 

\ahiu-Uh,Vh) + bhivh,p- Ph)\ 


sup 

Vh&Vh 




< Ch. 


This combined with sup < C'IIm - ^ Ch + C\\vh - Uh\\H^{nh) proves the claim. 

Consequently, we have 

( 4 . 6 ) Clip,, - 9 ft||L 2 (nh) < Ch-f C||i;,, -Uh\\m(Q.u) +Cph\\Viph - g,i)||L2(Q,^). 

Collecting the above estimates for /31,... ,135 and noting that ?7/i^||V(pft —(7 ,i)|||2(q^) = dhiph — qh^Ph — 
qh), we deduce 

( 4 . 7 ) I/3I < Ch^ + ^llr;,, - Uh\\]jna^) + ]^dhiph - qh,Ph - qh)- 
In the same way as we computed I3, one has 

h = Chiivh - u) - Uh, ph-^)+ Chiivh - u) - Uh, Xh - Ph) 


+ Ch{{u - Uh) ■ nh,Ph - A) 
+ Ch{{u - Uh) ■ Uh, Xh - Ph) + echiXh, Xh - Ph) 

— iChiph, Xh — Ph) ='- I41 + I42 + I43 + 744 -|- 745. 
By interpolation error estimates on F,,, we have 


|74i| < Ch\ I/ 42 I < Ch^/^WXh - Ph\W(T,^), 

I745I < Ce||A,, - piilUqrfc), 


whereas it follows from ( 4 . 2)3 that 


I 743 I < C/i^/^||{t — Mft||i/i(nh)) 


I744I < Ch\\Xh — phWh'^iTh)- 
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Noting that ||Ah — = Ch{\h — fJ-h, — IJ-h) and using Young’s inequality, we arrive at 

1/41 < Ch"^ + Ch^/e + Ch + Ce + Ch^/t 


+ "g ~ ’^hW'ni^Q,^) + — ^J'h, A/t — ^ih) 


(4.8) 


a, 


< C{h + e + h /e) + -^\\vh ~ UhWn^ 


(f^h 


4“ ^h: 


Combining (4.5), (4.7), and (4.8) with (4.4), we conclude the desired estimate for Ijiift — Uh\\H^{Qh)- 
The result for ||(Z/i — Ph\\L^{n^) follows from (4.6). This completes the proof of Theorem 


4.1 


□ 


Remark 4.2. According to the theorem, the best rate of convergence is obtained by choosing 

e = 0{h), which is not optimal. Let us highlight the reasons for this sub-optimality. First, as far as the 
variational principle is concerned, the most suitable regularity to work with for Xh would be 
instead of L‘^{Th) as presented above. However, it is not possible to extract this regularity from 144 
above (more precisely, Ij in the proof of Proposition 4.1) because ^ i7^/^(r^). Second, it is not 
trivial whether the following inf-sup condition would hold: 

Chivh ■ nh,Hh) 


(4.9) 


C\\f^h\\H~^/^(ru) < sup II ynh&^h- 

vheVn 


In the case fl/j = H, this condition is valid for a suitable choice of A/j. Qaglar and Liakos 1113 took 
advantage of this fact to derive the optimal rate of convergence 0 {h + e). 

5. Penalty FE scheme with reduced-order numerical integration 


In this section, we investigate problem (3.1) in which Ch is replaced with its reduced-order numerical 
integration c\ defined via the midpoint (barycenter) formula as follows: 

X! \S\H'ms)K'ms), X,^lGC{Th), 

S&Sh 

where IS”! denotes the area of S and ms is the midpoint of S when N = 2 (resp., the barycenter of S 
when N = 3). Because we exploit pointwise evaluation of functions, we assume higher regularity of the 
exact solutions as follows: 

uGW^’°°{n)^, pGW^’°°{n), AGtTi’°°(r), 

which implies / G L°°(H)^, g G ^ ^1,00 

Then the problem we propose reads: find (uh^Ph) G 14 x Qh such that 

f ah{uh,Vh) -\- bh{vh,Ph) + -ci{uh ■nh-g,Vh- nu) = {f,Vh)n^ + {f,Vh)rh 
(5.1) < e 

[ bh{uh,qh) = dh{Ph,qh) 

for all {vh,qh) G 14 x Qh- The well-posedness of this problem is obtained by the same manner as in 
Proposit ion |3.1[ We also find that its solution satisfies the following three-variable formulation as we 
derived (3.5): 

ah{uh,Vh) + bh{vh,Ph) +c\{vh ■ nh,Xh) = {f,Vh)n^ + {f,Vh)ru ^Vh G 14, 

bhiuh^qh) = dh{ph,qh) yqhGQh, 

cl(uh ■nh-g,fj-) = ecl(Xh, g) Vg G C{Th), 

where Xh is defined only on {ms '■ S G Sh} by Xh{ms) = \{uh ■ nt — 5)|ms- Likewise, the error analysis 
is mostly parallel to the arguments in Theorem |4.1| Thereby, in the sequel we only focus on what will 
change due to the replacement of Ch by c\. In doing so, first we observe that: 

Lemma 5.1. Let Vh G 14 and X G Then 

\ch{vh -nhA) - cl{vh • nh,X)\ < /i||T,i||Li(rol|A||i4.i,oo(n). 
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Proof. Since the midpoint (barycenter) formula is exact for affine functions, one obtains 

Chivh ■ nh, A) - cl{vh • n/i, A) = ^ Vh ■ UhiX - A(ms)) d'jh- 

JS 


SeSn' 


This combined with ||A — A(ms)||j;,cx>( 5 ') < ||A|j^^i,oo(Q)h (note that diam^ < h) concludes the desired 
result. □ 


Combining this lemma with the estimates of I 5 in the proof of Proposition 4.1 we see that (4.2)i 
remains the same even if we replace Ch by cj^. 


Next we consider the analysis of ( 4 . 2 ) 3 , namely, the estimates for / 7-/10 in the proof of Proposition 


4.1| To this end we introduce a semi-norm in Ah by 


By Lemma B.l and Cauchy-Schwarz inequality, /y is bounded by Ch^\^h\Ah if N = 2 and by Ch\^h\Ah 

we have \Is\ + j/gl < 


A.2 


if iV = 3. By the regularity assumption u,g G and by Proposition 

Ch^\f.ih\A,- We notice that Iio in the present situation is zero. Therefore, instead of (| 4 . 2[)3 we obtain 


(5.2) 


\cl{{u - Uh) ■ nh,Hh) + ecl{Xh, fj.h)\ < Ch^ IhhIah G Ah, 


where j = 2 if = 2 and j = 1 if = 3. 

Finally, we consider the estimates of 14 in the proof of Theorem |4.1[ This time we may choose 
gh = IhX as an interpolation of A. In view of the regularity assumption u G VF^’°°(fi) and A G 


and by virtue of Lemma B.l we have 


\hi\ < Ch^, \I42\ < Ch^lXh — gh\Ah^ Ihsl ^ Ch\\u — Uh\\H^(nh.)^ 

|.f45| < Ce\Xh — gh\Ah- 


By (5.2), J 44 < Ch^\Xh — fJ-hlAh- Consequently, instead of (4.8) we obtain 


\I 4 \ < C+ e + j e) + — llii/i — — fJ-h,Xh — gh)- 

From these observations, we arrive at the following result. 


Theorem 5.1. In addition to the hypotheses of Theorem \4. 1\ we assume that th e sol ution (u,p) of (2.5) 
possesses the VF^’°°(f2)^ x -regularity. Let {uh,Ph) be the solution of (5.1). Then there exists 

a constant C = C{N, LI, v, u,p), independent of h and e, such that 

(5.3) IIu — Mftll/yqnfc) + ll(p + fc/t) — PhllL^^nu) < C{h-\- ^), 

where j = 2 if N = 2 and j = 1 if N = 3. 

Remark 5.1. According to the theorem, choosing e = 0{h'^) gives us the optimal rate of convergence 
0{(h) when N = 2. When A^ = 3, at least we see that introduction of reduced-order numerical integration 
does not deteriorate the rate of convergence. Our numerical example given in the next section shows 
that it does improve the accuracy for N = 3 as, well. 


6. Numerical examples 


In the sequel, we refer to the schemes (3.1) and (5.1), i.e. without and with reduced-order numerical 
integration, as “non-reduced” and “reduced”, respectively. 
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6.1. Two-dimensional test. In this example, all the computations are done with the use of FreeFem++ 
m choosing the Pl/Pl element, i.e. I = 1, together with 77 = 0.01. Let be the unit disk, namely. 


^ = {(x.y) G K 

f,g,T given by 


< 1}. We consider the slip boundary value problem (1.1 1 for u = 1 and for 


/ = 


+7/2) -F 16y 
x{x‘^ + 

1 — x^ —xy 
-xy 1 - 


5 = 0, 

— 12xy 2 (a;^ — y"^] 
2 { x '^ — y^) —4a;y 


in which case we have the analytical solution u = {—y{x^ + y'^),x{x^ + y^))^, p = 8xy. We also 
introduce the solution of the no-slip boundary value problem denoted by (yno-siip^pno-siip)^ Namely, it is 
determined according to the same / as above and to the boundary condition = 0 on P. Figure 

| 6 . 1 | shows the velocity profiles of the two solutions; one notices the clear difference in their circulating 
directions and in the maximum modulus of velocity. 



Figure 6.1. Velocity profiles of u (left) and (right). 


On a mesh with h « 0.241, we computed numerical solutions of the slip boundary value problem 
using the non-reduced/reduced schemes for three choices of the penalty parameter e: 0 {h), 0{h?), and 
very small. The results are shown in Figure |6.2[ We find that the reduced scheme gives more robust 
and accurate approximate solutions. In fact, in case e = 0 {h), there are two spurious circulations inside 
n for the non-reduced scheme; we remark that refining a mesh and keeping e = 0 {h) did not suppress 
them. For smaller e the non-reduced scheme fails to capture the correct solution and seems to approach 
the no-slip boundary value problem. This behavior is somehow expected because letting e —>■ 0 in the 
penalty term of (3.1) implies (at least formally) the constraint Uh ■ rih = 0 on Th, which undesirably 


collapses into u/j = 0 on F^j as observed in Section However, the reduced scheme produces solutions 
which capture the slip boundary condition correctly for all e > 0 sufficiently small. Therefore, it is 


expected that the error bound (5.3) could be improved in such a way that the reciprocal of e would not 


appear (we conjecture that some inf-sup condition like (4.9) would be valid). 

Next we study the convergence property of the non-reduced and reduced schemes, whose solutions 
are denoted by (u^^,p^^) and respectively. We compare the convergence behavior of them 

with that of the numerical solutions computed with the Dirichlet boundary condition, which are denoted 
by (i-e., the boundary condition = u onTh is imposed). The linear solver is chosen as 

UMFPACK, and u — Uh and p — Ph are interpolated into the quadratic Hnite element space to compute 
errors in the associated norms. For convenience, we also report errors in the L^(f2;i)^-norm of velocity, 
and we remark that ||u||i 2 (f 2 ) « 0.886, ||u||//i(o) « 3.355, ||p||L 2 (a) « 2.894. The results are presented 
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Figure 6.2. Velocity profiles of numerical solutions computed with non-reduced (left 
column) and reduced (right column) schemes. For each row (top to bottom), e is chosen 
as O.l/i, 10“®, where h « 0.241. 
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Table 6.1. Convergence behavior of velocity in the H^{Uh)^-norm in the 2D test (top: 
e = O.l/i, bottom: e = O.l/i^). DOF means the number of degrees of freedom. 


h 

DOF 

||u- 

Rate 

Ik- Wftll 

Rate 

lk-<kl 

Rate 

0.316 

333 

1.043 

— 

0.575 

— 

0.464 

— 

0.165 

1182 

0.567 

0.94 

0.296 

1.09 

0.231 

1.07 

0.078 

4488 

0.310 

0.81 

0.145 

0.94 

0.114 

0.94 

0.045 

17391 

0.146 

1.37 

0.077 

1.30 

0.057 

1.28 

0.023 

69270 

0.074 

1.02 

0.039 

1.02 

0.028 

1.02 

0.012 

274956 

0.036 

1.16 

0.020 

1.05 

0.014 

1.10 

h 

DOF 


Rate 

Ik - <11 

Rate 

lk-<kl 

Rate 

0.316 

333 

1.683 

— 

0.479 

— 

0.464 

— 

0.165 

1182 

1.559 

0.12 

0.232 

1.12 

0.231 

1.07 

0.078 

4488 

1.618 

(<0) 

0.114 

0.95 

0.114 

0.94 

0.045 

17391 

1.412 

0.25 

0.057 

1.28 

0.057 

1.28 

0.023 

69270 

1.388 

0.03 

0.028 

1.02 

0.028 

1.02 

0.012 

274956 

1.293 

0.11 

0.014 

1.10 

0.014 

1.10 





h 


Figure 6.3. Convergence behavior of Iju — (left), ||u — (middle), 

\\p — Ph\\L^(nh) (right) in the 2D test. The triangles indicate the slope 0{h). 


6.1 


in Table 

and that \\u — W 


and Figure 6.3 for two choices of e. We find that ||u — = 0{h?) for e = 0{h?) 


) = 0(h) for • = NR, R and for e = 0(h), which cannot be explained by the 

e = 0 (h?) achieves the best accuracy is 


7i II (iih 

theory. Nevertheless, the fact that the reduced scheme with 
in accordance with the theoretical prediction. We see from Figure |6.3| that its accuracy in the energy 
norm is almost the same as that of 




6.2. Three-dimensional test. In this example, we focus on the verificati on o f convergence. Let D be 
the unit sphere, i.e., D = {(x, y, z) £ < 1}. We consider (1.1) for u = 1 and for /, g, r 

such that the analytical solution is 


/ 10x‘^yz{y - z)\ 

10y‘^zx(z — x) , p = 10xyz(x + y + z). 
\10z^xy(x - y)J 


u = 
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Table 6.2. Convergence behavior of velocity in the {Uh)^-norm in the 3D test (top: 
e = O.l/i, bottom: e = O.lh^). Itr means the number of iterations required for GMRES 
to converge. 


h 

DOF 


Rate 

Itr 


Rate 

Itr 


Rate 

Itr 

0.240 

1.11E-F4 

1.471 

— 

69 

1.048 

— 

70 

1.268 

— 

76 

0.113 

1.12E-F5 

0.656 

1.08 

238 

0.638 

1.06 

279 

0.574 

1.06 

349 

0.075 

3.42E-F5 

0.454 

0.89 

352 

0.448 

0.86 

352 

0.405 

0.85 

393 

0.062 

6.59E-F5 

0.372 

1.08 

469 

0.367 

1.07 

473 

0.326 

1.16 

751 

0.052 

1.17E-F6 

0.306 

1.05 

655 

0.303 

1.04 

658 

0.266 

1.10 

790 

0.045 

1.88E-F6 

0.263 

1.04 

901 

0.261 

1.03 

899 

0.227 

1.09 

1979 

0.039 

2.39E-F6 

0.238 

0.87 

1749 

0.236 

0.87 

1638 

0.205 

0.88 

5179 

h 

DOF 


Rate 

Itr 

Ih-Mftll 

Rate 

Itr 


Rate 

Itr 

0.240 

1.11E-F4 

1.590 

— 

77 

1.350 

— 

81 

1.268 

— 

76 

0.113 

1.12E-F5 

0.877 

0.79 

270 

0.579 

1.13 

304 

0.574 

1.06 

349 

0.075 

3.42E-F5 

0.630 

0.80 

467 

0.405 

0.87 

646 

0.405 

0.85 

393 

0.062 

6.59E-F5 

0.575 

0.49 

742 

0.327 

1.15 

782 

0.326 

1.16 

751 

0.052 

1.17E-F6 

0.534 

0.41 

1111 

0.271 

1.02 

1348 

0.266 

1.10 

790 

0.045 

1.88E-F6 

0.493 

0.54 

1735 

0.231 

1.08 

2175 

0.227 

1.09 

1979 

0.039 

2.39E-F6 

0.477 

0.29 

2103 

0.201 

0.89 

2600 

0.205 

0.88 

5179 


We remark that ||u||i 2 (Q) « 0.708, ||M||_Hi(a) « 4.943, ||p||L 2 (n) ~ 1.043. The computations are done 
with the use of FEniCS |23] (combined with Gmsh [12] for obtaining meshes) choosing the PI/PI element 
with r] = 0.1. The linear solver is GMRES, preconditioned by incomplete LU factorization, with the 
restart number 200 and with the relative tolerance 10“®. As in the previous example, u — ut and p — pt 
are interpolated into the quadratic finite element space to compute their norms. The results are reported 

Except case of the non-reduced scheme with e = O(h^), the errors seem to 


6.2 and Figure 6.4 


in Table ^^ ^^ 

converge at the rate 0{h), which is better than the theoretically predicted one Our opinion 

is that the interior errors would be dominant with the resolution of meshes considered here and that 
the suboptimal rate 0{h^^^) would be observed only for a very fine mesh. However, from the results 
we infer that the reduced-order numerical integration is also effective for the case = 3, in which the 
choice of e = O^h?) may be recommended (although this was not justified by a rigorous proof). With 
this choice, as in the two-dimensional case, the accuracy in the energy norm is comparable with that of 


6.3. Affect of penalty parameter on linear solvers. It is known that the use of too strong penalty 
(i.e. small e in our case) would lead to ill-conditioned problems (see e.g. 0), which could deteriorate the 
performance of linear solvers. Therefore, we examine a condition number of the matrix A obtained from 
our penalty FE scheme, in particular, its dependency on the penalty parameter e. For this purpose, 
we fix the mesh h « 0.113 in the 3D test, varying e from 100 to 10“®. Moreover, we consider only the 
reduced scheme since the behavior was similar for the non-reduced one. The condition number is then 
estimated by the Matlab function condest (A). We also report the number of iterations required for 
GMRES and BiGGSTAB to c onve rge, which were preconditioned by incomplete LU factorization. The 
results are presented in Table 


6.3 


It seems that the condition number grows at the rate 0(e“^), which 
is faster than the one explained in [^ p. 532]. One also notices that GMRES and BiGGSTAB failed to 
converge when e < 10“®. We remark that, even for such small e, sparse direct solvers like UMFPACK or 
MUMPS were able to solve the linear system (apparently with no problem). Although we do not have 
a good explanation of this phenomenon, it tells us that one should carefully choose e in order to assure 
both the accuracy and the numerical stability, especially when one wants to invoke iterative methods 
for solving linear systems obtained from the penalty method. 
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10 *^ 


10 ' 

o 

LU 


10 '‘ 


10 -' 


Reduced(eps=0.1h2) - 
Reduced(eps=0.1ho - 
Non-Reduced(eps=0.1 hp) - 
Non-Reduced(eps=0.1 h ) - 
Dirichlet 



10“ 

o 

LU 


Reduced(eps=0.1hi) ■ 
Reduced{eps=0.1hJ ■ 
Non-Reduced{eps=0.1 hp) - 
Non-Reduced{eps=0.1 h ) ■ 
Dirichlet 



10‘ 

h 


10' 


10 ‘ 

h 



Figure 6.4. Convergence behavior of |ju — Wft||z, 2 (nh) (left), \\u — Uh\\H^{nh) (middle), 
\\p — Ph\\L^{nh) (right) in the 3D test. The triangles indicate the slope 0{h). 


Table 6.3. Penalty parameter, condition number, and number of iterations for GM- 
RES (the 4th and 5th columns) or BiCGSTAB (the last column) to converge in the 3D 
test with the mesh h » 0.113 (DOF is 112476). The absolute and relative tolerances 
are set to 10 “^° and 10 “®, respectively. 


e 

condest(A) 

Rate 

Itr (restart=30) 

Itr (restart=200) 

Itr (BiGGSTAB) 

1.0E-P2 

2.36E-P6 

— 

655 

182 


1373 

l.OE-Pl 

2.27E-P6 

(< 0 ) 

672 

191 


165 

l.OE-PO 

2.45E-P6 

0.03 

733 

195 


516 

l.OE-1 

3.53E-P6 

0.16 

392 

165 


264 

l.OE-2 

1.64E-P7 

0.67 

480 

139 


152 

l.OE-3 

1.46E-P8 

0.95 

539 

195 


894 

l.OE-4 

1.56E-P9 

1.03 

1432 

352 


2888 

l.OE-5 

1.54E-tll 

2.00 

28162 

370 


2293 

l.OE -6 

1.54E-tl3 

2.00 

(not converged) 

(not converged) 

(not 

converged) 

l.OE-7 

1.54E-tl5 

2.00 

(not converged) 

(not converged) 

(not 

converged) 

l.OE -8 

1.54E-tl7 

2.00 

(not converged) 

(not converged) 

(not 

converged) 


7. Conclusion 

We investigated the Pl/Pl and Plb/Pl finite element approximations for the Stokes equations subject 
to the slip boundary condition in a domain with a smooth boundary. The constraint u ■ n = g on T 
is relaxed by using the penalty method, which enables us to avoid a variational crime and makes the 
numerical implementation easier. We developed a framework to address the difficulty due to D 7 ^ D/j 
and successfully applied it to establish error estimates of the finite element approximation. The use 
of reduced-order numerical integration together with e = 0{h?) in the penalty term improves the 
accuracy, which was theoretically justified for N = 2 and was numerically confirmed for TV = 3. In fact, 
we observed that the accuracy in the energy norm was comparable with that of numerical solutions 
subject to Dirichlet boundary conditions. 
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Appendix A. Transformation between r and F/j and related estimates 

Let ri be a bounded domain in Its boundary F is assumed to be C^’^-smooth, namely, there 
exist a system of local coordinates {{Ur,yr,^r)}^i and positive numbers a, (3 such that: 1) {Ur}^i 
forms an open covering of F; 2) i/r = (Vri, ■ ■ ■ ,yrN-i,yrN) = [y'r^yrN) is a rotated coordinate of the 
original one x, that is, yr = ArX for some orthogonal transformation 3) ipr S gives a graph 

representation of F n Ur, where A^. := {y'^ e : \y'^\ < a}, that is, 

F n t/r = {i/r e : y'r^ A^ and yrN = ^riv'r)}, 

nnUr = {yr : y'r€ Ar and Priy'r) < yrN < Priy'r) + /3}, 

n'' nUr = [yr e : y'r & Ar and (friy'r) - P < VrN < Pr{y'r)}- 

Since F is compact, there exists ho > 0 such that for any x G F the open ball B{x',ho) is contained 
in some local coordinate neighborhood Ur- According to the fact that C^’^(Aj.) = IF^’°°(Aj.), the 
derivatives of pr are bounded up to second order, i.e., 

||</5r|U“(A,.) < Co, ||VVr|U“(A,.) < Cl, || 1| (A,.) < C2, 

where Co, Ci, C 2 are constants independent of r, and V' means Vy'. 

The smoothness of F is connected with that of the signed distance function d{x) defined by 

—dist(x, F) if X G fl, 
dist(x, F) if X G 

We collect several known properties on d{x) below. For the details, see e.g. [131 Section 14.6] or [H 
Section 7.8]. Let F(^) := {x G R'^ : ld(x)j < i5} be a tubular neighborhood of F with width 25. Then 
there exists <5 depending only on the curvature of F such that for arbitrary x G F((5) the decomposition 

(A.l) X = 7r(x) + (i(x) n(7r(x)), 7r(x) G F, 

is uniquely determined. Here, n is the outer unit normal held dehned on F, which coincides with Vdjr- 
We extend n from F to F(5) by n(x) = Vd(x), which also agrees with n(7r(x)). The fact that F is 
C^’^-smooth implies that d G (7^’^(F((5)), n G C'°’^(F((5)), and tt G C°’^(F((5)). We call tt : F((5) —>■ F the 
orthogonal projection onto F, in view of its geometrical meaning. We may assume that 

IMIlL“(r(5)) < Co , ||V(iljLoo(r(5)) < Ci , j|V^dljLoo(r(5)) < C2 , 

where we re-choose the constants Co, Ci, C 2 if necessary. 

Now we introduce a regular family of triangulations {Th\h\.o of H in the sense of Section 
we denote by Sh the boundary mesh inherited from Tn, and we set H/j = UreTfeT’ and F/j = \Js^Sh^- 
In order for F/j to be compatible with the local-coordinate system {{Ur,yr,Pr)}^=i, we assume the 
following: 

1) the mesh size h is less than min{/io, 1}; 

2) for every r = 1,..., M, F/i n 17^ is represented by a graph {(?/', prh{y'r)) G R'^ : y'r G A^}; 

3) every vertex of S' G 5?^ lies on F. 

From these we see that every S G S/i is contained in some Ur and that prh is a piecewise linear 
interpolation of pr- As a result of interpolation error estimates, for arbitrary r we obtain 

ll7’r?i||L°°(A,.) S Co, \\pr — Prh\\L°°{Ar) S Cosh^, 

II VVrIi||L°°(A^) < Cl , ||V'((pr — Prh )\\ L ^( Nr ) < 

where we re-choose Co, Ci if necessary and the subscript E refers to “error”. In the following, h is made 
small enough to satisfy 2CoEh^ < min{ho, ^}, which in particular ensures that tt is well-defined on F^. 
We may assume further that 7r(5') is contained in the same Ur that contains S. 

Based on the observations above, we see that the orthogonal projection tt maps F/j into F. It indeed 
gives a homeomorphism between F?i and F, and is element-wisely a diffeomorphism, as shown below. The 
representation of a function /(x) in each local coordinate (Ur,yr,Pr) is defined as f{yr) ■= fiA~^yr). 


3.1 


As before, 
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However, with some abuse of notation, we denote it simply by f{yr)- Then, the local representation of 
the outer unit normal n associated to T is given by 


n{yr,<fr{yr)) = 


Kr{y'r) 


VVr(2/;) 


yr ^ 


where Kriy'r) '■= \/l + |VVr-(2/r)P- We are ready to state the following. 

Proposition A.l. If h > 0 is sufficiently small, then Trlr^ is a homeomorphism between T^ and T. 

Proof. Let us construct an inverse map tt* : T —>■ T^t of ttIt^ which is continuous. To this end, we fix 
arbitrary a: € T and choose a local coordinate (Ur,yr,<Pr) of x such that B{x;ho) C Ur- For simplicity, 
we omit the subscript r in the following. In view of the definition of 7r(a;) by (A.l), we introduce a 
segment given by 


+ 


vv(2/'; 

-1 


|t| < 2CoEh^- 


^<p{y')J ^ K{y') 

To each point on this segment we associate its height H{f) with respect to the graph of tph^ that is, 

m = viy') + VV(»')) . 

Then we assert that < 0 and that iL(—2 Coe/i^) > 0, HiflC^Eh^) < 0. 

To prove the first assertion, letting 

we have ^H{t) = + VV(2/0 ' One sees that 

V'ph{Y') = vV(?/') - VV(j/') + VV(W) - VV(W) + V'^h{Y') 

=-.V'v{y') + hPl2, 

where Ii and I 2 satisfy 

|/i| < ||V'VIU-(A)|W - y'\ < C2 ■ 2 CoEh\ I/2I < CiEh. 

Then it follows that 

1 + VV(y') • VV/.(W) > K{y'f - Ci{2C2CoEh^ + Oi^h). 

From this we have ^H{t) < 0 provided Ci{2C2CoEh^ + CiEh) < 1/2. For the second assertion, one 
finds that 

H{-2CoEh^) = + Piy') - Ph{Y') 

= - <f{Y') + ip{Y') - iphiY'), 

where Y' is y' — By Taylor’s theorem, there exists some 6 e (0,1) such that 

^(W) - = VV(j/') • (W - y') + l(Y'- y'fVWy,+eiY'-y'){Y' - y') 

- + 2C2ClEh\ 

By the definition of K{y) and by \<p{Y') — iph{Y')\ < Cob/i^, we obtain 

H{-2CoEh^) > 2K{y')CoEh^ - ‘2C2C^Eh^ - C^Eh^ 

> 2CoEh^ - ‘2C2ClEh^ - CoEh^ = CoEh^il - 2C2CoEh^), 

which implies that iL(—2 C'o_e/i^) > 0 provided C 2 C'oEh^ < 1/4. In the same way, the last assertion 
iL(2C'o£)h^) < 0 can be proved. 
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From these assertions we deduce that there exists a unique t*(x) £ [—2Cob^^, 2Cob^^] such that 
H{t*{x)) = 0. Consequently, the map tt* rF —J-F^; x ^ x + t*{x)n{x) is well-defined. A direct 
computation combined with the uniqueness of the decomposition (A.l) shows that tt* is the inverse of 
7r|r,j. The continuity of tt*, especially that of C, follows from an argument similar to the proof of the 
implicit function theorem (see e.g. [20l Theorem 3.2.1]). □ 


Proposition |A.1| enables us to define an exact triangulation of F by 

7T{Sh) = {^(5) : 5 e 54. 

In particular, we can subdivide F into disjoint sets as F = UseSh Furthermore, for each S £ Sh we 

see that S and 7r(S') admit the same domain of parametrization, which is important in the subsequent 
analysis. To describe this fact, we choose a local coordinate {Ur, yr, ^r) such that Ur D S LI ^{3), and 
introduce the projection to the base set br : —>■ by br{yr) = y'r- The domain of parametrization 

is then defined to be S' = br{TT{S)). We observe that the mappings 

$ : S" Tr{S); 4 (4, ipriy'r))'^, 

^h--S' ^ S; y'r ■n*{y'r,y}r{y'r)) = ^(4) + «(^(2/r)). 

are bijective and that $ is smooth on S'. If in addition especially t*, is also smooth on S', then $ 
and $ may be employed as smooth parametrizations for S and Tr{S) respectively. The next proposition 
verifies that this is indeed the case. 


Proposition A. 2. Under the setting above, we have 

PIIl“(S') < CoBh'^ , |iV't||i=o(5/) < CiBh, 

where Cqe and Gib are constants depending only on N and F. 


Proof. Since the first relation is already obtained in Proposition A.l with Cqb = 2Coe, we focus on 


proving the second one. For notational simplicity, we omit the subscript r and also use the abbreviation 
di = ^ {i = 1) • • ■ j N). The fact that t* is differentiable with respect to y' can be shown in a way 
similar to the proof of the implicit function theorem. Thereby it remains to evaluate the supremum 
norm of V't in S', which we address in the following. 

Recall that t*{y') is determined according to the equation 

(A.2) i{y') = (p{y') - cph{y' + i{y')V(p{y')), 

where we have set i{y') := t* {y')/K{y'). Because V't* = KV'i + ^ v f* ^ follows that 

l|vr< (1 + Ci)||V't|Uo=(sq + C2C^CoEhG 

Therefore, it suffices to prove that ||VT||icxj( 5 /) < Ch] here and hereafter C denotes various constants 
which depends only on N and F. 

Applying V' to (A.2) gives 

(A.3) (1 + VV(2/') • VV4h^'))V'f = VV(2/') - Viph{Y') + i{y')V^ip{y')Vy>h{Y'), 

where Y' := y' + i{y')\7'(p{y'). By the same way as we estimated Ii and I 2 in the proof of Proposition 
ED we obtain 

|VV(j/') - VV4W)I < C2CoEh^ + CiEh < Ch, 


Also we see that 


1 + VV(y') • VV4W) > K{y'f - Ci{C 2 CoeIY + CiEh) > -. 


\i{y')V^<p{y')VMY')\ < CoEh^ ■ C 2 C, < ChG 


Combining these observations with (A.3), we deduce the desired estimate ||VT||ioo(s/) < Ch. 


□ 


Remark A.l. Let F £ Since iph is linear on h{S), further differentiation of (A.3| gives us 

t* £ C^’^{Tr{S)); in fact we have ||< G 2 e- This implies that 7r|s is a C^’^-diffeomorphism 
between S and 7r(S'). However, since V(fh is smooth only within b{S), tt is not globally a diffeomorphism. 
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Now we give an error estimate for surface integrals on T and F/j. Heuristically speaking, the result 
reads Idj — d'jhl < 0{h?), which may be found in the literature (see e.g. m)- Here and hereafter, we 
denote the surface elements of F and F/i by dy and d^h , respectively. 


Theorem A.l. Let S € Sh and / be an integrable function on S. Then we have 


'<s) 


fd'y- / foTTd'jh 


< Ch^ 


\f\dl, 


where C is a constant depending only on N and F. 

Proof. Let {Ur,yr,Tr) be a local coordinate that contains S U We omit the subscript r and use 

a 

dyi 


the abbreviation di = {i = 1,..., N). We represent the surface integral using the parametrization 


$ as follows: 


T(S) 


fd'y= [ f{<^{y'))VdetGdy', 
Js' 


where G = (Gij)i<ij<Ar-i denotes the Riemannian metric tensor given by Gij = 9^$ • (dot means 
the inner product in M^). Similarly, noting that tt o ih/j = $, one obtains 


[ foTrd-fh = [ f{<^{y'))\/det Gh dy', 

JS JS' 

where Gh is given by Gh,ij = di^h ■ dj^h- Then we assert that: 

\\Gh-G\\L^^S')<Ch^- 

To prove this, noting that = $ + t*n o $, we compute each component of Gh — G as follows: 


Ghhj - Gij =di^- dji^h - ^') + • di{^h - $) + di{^h - $) • dj{^h - $) 

= di^ ■ dj{t*n o $) + dj^ ■ di(t*n o $) + di{t*n o $) • dj(t*n o $) 
=: Ii +12 + 13- 


For Ii, we notice that 9^$ is a tangent vector so that 9^$ • n o <i) = 0. This yields 


h = 9i$ • t*dj{n o $) = t*di^ ■ (djn + djipd]s[n)\^, 


which is estimated by C'o_e/i^( 1 + C'i)(C'2 + G 1 G 2 ) thanks to Proposition 
same manner. To estimate /a, we observe that 


A.2 


I 2 can be bounded in the 


di[t*n o $) = {dit*)n o $ + t* {pin + diipd^n) |$, 

which is bounded by CiEh + C'o_e/i^(C' 2 + C'iC'2) < Gh. Similarly one gets \dj{t*n o $)| < Gh, hence it 
follows that I/3I < Ch"^. Therefore, \Gh,ij — Gp < CK^, which proves the assertion. 

Now we use the following crude estimate for perturbation of determinants (cf. [171 equation (3.13)]): 
if A and B are N x N matrices such that \Aij\ < a and \Bij\ < b for all i,j, then 

jdet (A + B)- det A| < N\N{a + b)^-^b. 


Combining this with the assertion above and also with ^/a — y/j3 = {a — (3)/{^/a + yp), we obtain 

WVd^-PdetGhh^iS') < Ch^. 

In addition, note that VdetG = a/ 1 + jVVP > 1- Consequently, 


t(S) 


fd-y- / foTTdjh 


<Gh^ f \fmP))\Vd^dy'= Ch^ f l/jdy 

JS' Jtv{S) 


which proves the theorem. 


□ 














22 


TAKAHITO KASHIWABARA, ISSEI OIKAWA, AND GUANYU ZHOU 


Remark A.2. Adding up the results of the theorem for all S G Sh yields 


fd'y- foTTdjh 


<Ch^ J^\f\dj. 


It also follows that | /p / o nd^hl ^ C Jj. \ f\ d'j. Choosing in particular \f\P as the integrand gives 
ll/o7r||ip(r^) < C'i/P||/||LP(r) for p G [l,oo]. 

Let / be a smooth function given on F. Then its transformation to F^ is defined by / o tt. However, 
if / is extended to a neighborhood of F, e.g. to F((5), then we may also consider /’s natural trace on F/j. 
The next theorem provides error estimation of these two quantities. 

Theorem A.2. Let f G VF^’^(F(Ji)), where p G [l,oo] and G 2C'o£;/i^]. Then, 

\\f - f° T^\\Lp{rh) < C5\ ^^^|!/|| wi.p(r(5i)), 
where C is a constant depending only on p, N, LI. 

Proof. Since F(5i) = UseSh7r(5',5i) = {a; G F(5i) : 7r(a;) G tt{S)} denotes a tubular 
neighborhood of S, it suffices to prove that 


(A.4) 


[\f-fo7 T\Pd-f<C5P-^[ \Vf\Pdx \/S 
Js 


€ Sh- 


To this end, using the notation in Theorem A.l we estimate the left-hand side of (A.4) by 


[ I/-/oTrl^’dy = [ 1/ o / o $|p det Gh dy 
Js J S' 

<C [ \fo<i>h-fo<i>\Pdy'. 

Js’ 


Here, for fixed y' G S' we have 


fi'^hiy')) - /(4>(y')) = J + ^i‘^h{y') - $(?/'))) ds 

= ® i*{y')n{Hy'))) ds 


= / t*{y)n{‘^{y))-L/f{<^(y) + st*{y)n{<^{y)))ds 

JO 


u($(i/')) • V/($(y') -f tn{<^{y'))) dt. 


Because \t*{y')\ < CoEh^ < <5i, it follows that 

\f{<^h{y'))-fmy'))\< r \Vfmy')+tnmy')))\dt 

J-Si 


< (25i 


A-l/p 


/-5i 


\ 1/P 

\Vfmy')+tnmy')))\^dt\ , 


where we have used Holder’s inequality. Consequently, 

(A.5) [ \f - f o Tr\P d'j < CSl~^ [ \Vf{^{y')+tn{^{y')))\^dy'dt. 

Js dS'x[-<5i,5i] 

On the other hand, we observe that the A-dimensional transformation 

(A.6) 4': 5" X 7r(5', ^i); {y',t) i-G ^{y') + tn{^{y')) 

is bijective and smooth. Application of this transformation to the right-hand side of (A.4) leads to 

f \L7f\Pdx= f \Wfmy')+tnmy')))f\detJ\dy'dt, 

J-rviS.Si) dS'x[-5i,5i] 
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where J = (9i$ + t9i(no$), • • • ,dN-i^ + tdN-i{no^),no^) denotes the Jacobi matrix of ih. Letting 
,dN- 1 $, n o $), we find that 

II J - J||l=o(S'x[-<5i,5i]) < CSi, 

because \tdi{n o $)| = \t{din + di(pdNn)\,;,\ < Si{C 2 + C'iC' 2 ). This implies 

||det J-det J||l=o(5.,x[_5i,5i]) < C 61 , 

which combined with det J = K{y') > 1 yields det J > 1/2 if h is sufficiently small. Therefore, 

(A.7) [ \Vf\Pdx>l [ \Vfmy')+tnmy')))\^dy'dt. 

Jtt{S,Si) jS' x[-Si,Si] 

The desired estimate (A.4) is now a consequence of (A.51 and (A.7). This completes the proof. □ 

Finally, we show that the L^-norm in a tubular neighborhood can be bounded in terms of its width. 
Such estimate is stated e.g. in |35l Lemma 2.1] or in [28l equation (3.6)]. However, since we could not 
find a full proof of this fact (especially for = 3) in the literature, we present it here. 


Theorem A.3. Under the same assumptions as in Theorem A.2, we have 

ll/llLp(r(5i)) < C'(^i||V/||LP(r(5i)) + ^i^^ll/llLp(r)), 
where C is a constant depending only on p, N, fl. 

Proof. We adopt the same notation as in the proofs of Theorems |A.1| and ] A. 2[ Then it suffices to prove 
that 


f \f\^dy<clslf \Vf\Pdy + 6,f \f\Pdj] WS 
J7r{SJi) \ JTr{S,Si) J-KiS) J 


€ Sh- 


To this end, using the transformation di given in (A.6) we express the left-hand side as 
\f\Pdy= [ \f{'i){y',t))\P\detJ{y',t)\dy'dt 


r{S,Si) 


IS'x[-S,.S,] 


<C 

Js'x[-(5i.5i] 
=:/i+/2. 


For Ii, we see from the same argument as before that 


iifinv', t)) - fww + umy'w) dy'dt 


\f{d>{y'U))-fWW<{25,) 


p-i 






\Vfmy') + snmy')))fds, 


which yields 


\Ii\<CSl[ \Vfid>iy\s)f\detJiy\s)\dy'ds = C6l[ iVffdy. 

Js'x[- 5 i, 5 i] JT^iS.Si) 

For I 2 , it follows that 

\I2\ = 2CS, [ IfmWdy' <2C6i [ \fmy'))\PV^dy' = 2CS, [ l/^dq. 
JS' Js' Jtv{S) 

We have thus obtained the desired estimate, which completes the proof. 


□ 
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Appendix B. Error of n and Uh 

Let us prove that \n o tt — nfi\ < 0{h) on and also that, when = 2, it is improved to 0{h?) if 
the consideration is restricted to the midpoint of edges. 


Lemma B.l. Let n and Uh be the outer unit normals to T and respectively. Then there holds 

(B.l) ||R0 7r-n;,||i=°(rfe) < C'/i- 

If in addition N = 2, T ^ and ms denotes the midpoint of S G Sh, then 

(B.2) sup \n o 7r(ms) — nh{ms)\ < Ch^. 

seSh 

Here, C is a constant depending only on N and F. 


Proof. Let S G Sh he arbitrary and let {Ur,yr,(fr) be a local coordinate that contains S U 7r(S'). We 
omit the subscript r in the following. One sees that n and Uh are represented as 


n{y',T{y')) = 


1 


vv' 


nh{y',Th{y')) = 


1 


^/T+WW 


vV/t 


\/i + |vvP 

A direct computation gives 

(B.3) \n{y\ Lp{y')) - nh{y\ ‘Phiy'))\ < 2|V'(v?(y') - iph{y'))\ < 2CiEh. 

This combined with the observation that 

\noTT{y',Lph(:y')) - n{y',ip{y'))\ = \no TT{y',iph{y')) - no Trijj',ip{y'))\ 

< Ci\TT{y',(ph{y)) - 7r(?/',(p(y'))| 

< C'il|V7r||ioo(r(5))|(/?ft,(y') - ‘p{y')\ 

< OillV7r||2,cx=(r(5))C'oEft.^ 


2/'e A. 


proves (B.l). When N = 2 and F G by using Taylor expansion, we fin d th at ms is a point of 

super-convergence such that \ ‘^{h{ms)) — ^{h{ms))\ < Ch'^. This improves (B.3| to 0{h?), and thus 
(B.21 is proved. □ 


Remark B.l. In the case A^ = 3, it is known that the barycenter of a triangle is not a point of super¬ 
convergence for the derivative of linear interpolations; see 
only for A^ = 2. 


p. 1930]. For this reason, (B.2) holds 
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